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o Abstract 



11 We use the non-linear reduced- MHD code JOREK to study ELMs in the geometry of the AS- 

12 DEX Upgrade tokamak. Toroidal mode numbers, poloidal filament sizes, and radial propagation 

13 speeds of filaments into the scrape-off layer are in good agreement with observations for type-I 
u ELMs in ASDEX Upgrade. The observed instabilities exhibit a toroidal and poloidal localiza- 

15 tion of perturbations which is compatible with the "solitary magnetic perturbations" recently dis- 

16 covered in ASDEX Upgrade [R.Wenninger et.al., Solitary Magnetic Perturbations at the ELM 

17 Onset, Nucl.Fusion, accepted, preprint at http://arxiv.org/abs/1202.3 603]. This 
is localization can only be described in numerical simulations with high toroidal resolution. 
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ig 1. Introduction 



20 Edge localized modes (ELMs) accompany the high-confinement mode (H-mode) in tokamak 

21 fusion plasmas. As they could cause a potentially destructive heat-load on divertor plates and 

22 wall structures in large fusion devices [1, 2], detailed understanding of these instabilities and 

23 of mitigation-techniques is important for a successful operation of ITER. Non-linear MHD- 

24 simulations with the JOREK code [3-9], which is also used for the present article, and other 

25 codes like BOUT++ [10], NIMROD [11], or M3D [12] can make an important contribution 

26 after successful benchmarks with measurements in existing tokamaks. 

27 In the present article, a comparison between simulations with the non-linear finite-element code 

28 JOREK [13] and observations in the ASDEX Upgrade tokamak [14] is started. We concentrate 

29 on the early phase of ELMs. JOREK solves the reduced MHD equations in realistic X-point 

30 geometries as described in Section 2. ASDEX Upgrade is equipped with a unique set of edge di- 

31 agnostics that allows to investigate ELM crashes with high spatial and temporal resolutions [15]. 

32 This provides excellent possibilities for theory-experiment comparisons. Emphasis is put on 

33 simulations with high toroidal resolution (many toroidal modes at toroidal periodicity 1) to treat 

34 the coupling between various toroidal modes properly. This way, aspects can be identified that 

35 are described well already at low toroidal resolution (few toroidal modes at a high toroidal peri- 

36 odicity) while others are influenced significantly by the non-linear toroidal mode-coupling. 

37 The article is structured as follows. Section 2 describes the non-linear MHD-code JOREK. 

38 Physical parameters and technical details of the numerical simulations are given in Section 3. 

39 Our observations and findings made in the simulated instabilities are presented in Section 4. 

40 Subsequently, Section 5 describes how these results compare to experimental measurements. 

41 Finally, Section 6 summarizes and gives a brief outlook. 



2. JOREK Code 



43 The simulations are carried out with the single-fluid reduced-MHD model of the JOREK code. 

44 Section 2. 1 describes the equations solved in this model. For more details on the derivation, refer 

45 to Reference [16] and Appendix A. Spatial and temporal discretizations are briefly addressed in 

46 Section 2.2. 
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Table 1 : The normalization of quantities in JOREK is listed, it corresponds to choosing scale 
factors Bo = 1 T and Rq = 1 m. Variable names with subscript "SI" denote quantities in 
SI units, while variables without this subscript are the ones used in JOREK. In the pre- 
sented simulations, no = 6 • 10 19 m~ 3 and po = 2 • 10~ 7 kgm -3 . The magnetic constant 
is denoted and the Boltzmann constant kg. 
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47 2.1. Reduced-MHD Equations 



48 Seven physical variables are treated: poloidal flux stream function u, toroidal current density 

49 j, toroidal vorticity CO, density p, temperature T, and velocity vii along magnetic field lines. The 
so normalization of the relevant quantities is listed in Table 1. 

51 Variables j and co are connected to *F and u by the definition equations 



;=A«*=^,. (R -v^)= R A(if) + g, (i, 



<0 = v ^"-r1r{ r 7r) + ^ (2) 



Id/ du\ d u 



52 where V po i denotes the del-operator in the poloidal plane, 7? the major radius, and Z the vertical 

53 coordinate. The time-evolution of the remaining five free variables is described by the following 

54 set of equations (called physics-model "302" in JOREK): 



d^> r . du 

^ = Vj -R[u^-F -, (3) 



dp 

dt 



-V-(pv) + V-(D ± V ± p)+5 p , (4) 



dT 

p-^ = -p\-VT-(K-\) P V-\ + V- (K ± V ± T + K ll V ll T)+S T , (5) 
e r Vx Jp|^ = -p(v-V)v-Vp+jxB + vAv|, (6) 

B-|p^ = -p(vV)v-V/j+jxB + vAv|. (7) 

55 In every time-step, Equations (1-7) are solved simultaneously in weak form as a large sparse 

56 implicit system of equations. The velocity vector is defined as 



-RVuxeA+vu B, (8) 
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Table 2: The toroidal mode numbers resolved in simulations with different periodicities are 
listed. 



Periodicity 


Resolved n-modes 


8 


0,8,16 


4 


0,4,8,12,16 


2 


0,2,4,..., 12, 14, 16 


1 


0,1, 2,. ..,14,15, 16 



57 the magnetic field vector as 



B = (F ^ + V^x^)/R, (9) 

58 the pressure is p = pT , and K = 5/3 denotes the ratio of specific heats. Here, denotes the 

59 normalized toroidal basis vector. The toroidal magnetic field = Fo /R is fixed and cannot 
eo change with time. The poloidal velocity, i.e., the velocity vector in the poloidal plane, is denoted 

61 Vpoi =. The parallel gradient is given by V| | = b(b-V), where b = B/|B|, and the perpendicular 

62 gradient by Vj_ = V - V 1 1 . The Poisson bracket [u , *F] is defined as ff f| - J§ f| . Note, that the 

63 poloidal components of the velocity in this set of equations are determined only by the ExB- 

64 drift term. As a result, u acts as a velocity stream function and (except for a factor Fo) also as 

65 electric potential. 

ee Ideal-wall boundary conditions are implemented where the boundary of the computational do- 

67 main is parallel to the magnetic flux surfaces. At the divertor targets, where the flux surfaces 

ea intersect the computational boundary, modified Bohm boundary conditions apply [4, 17]. 



eg 2.2. Discretization 

70 The poloidal plane is discretized by 2D Bezier finite elements with four degrees of freedom per 

71 grid node and physical variable [13], while a Fourier decomposition is applied toroidally. The 

72 number of toroidal modes resolved in the simulations and the assumed toroidal periodicity of the 

73 system can be chosen separately. A periodicity equal to one means that the solution is computed 

74 for the whole torus. For larger periodicities, only a toroidal section of the torus is resolved. The 

75 modes included in the presented simulations are listed in Table 2. 

76 The temporal discretization is performed by a fully implicit second-order linearized Crank- 

77 Nicholson scheme [18]. In the resulting large sparse system, all physical equations and all 



5 



Table 3: Core values for plasma resistivity and viscosity are listed for the simulations denoted 
eta5 and eta6. Both quantities are modeled with a r~ 3 / 2 dependence and are chosen 
significantly larger than in experiments due to computational restrictions. In ASDEX 
Upgrade, the core resistivity is typically about 10~ 8 £lm. 



Run 


T]si [flm] 


Vsi [ m 2 /s] 


eta5 


5 x lfr 5 


7.5 x 1(T 5 


eta6 


5 x 1(T 6 


7.5 x 10~ 6 



78 toroidal harmonics are coupled. It is solved by an iterative GMRES -method, where a physics- 

79 based preconditioning is applied at the beginning of each GMRES solver step. In the precondi- 
ao tioning, the coupling between the sub-matrices corresponding to individual toroidal harmonics 
si is neglected which allows to solve each sub-system separately. This is performed using the direct 
82 solver PaStiX [19]. 



83 3. Simulations 



84 Simulations of edge-localized modes are one of the most challenging tasks in fusion MHD 

as numerics. The problem must be treated in realistic X-point geometry as the mode-affected 

ae region extends from inside the H-mode pedestal into the scrape-off layer, the vacuum region, 

87 and to the divertor legs. High spatial resolutions in all dimensions are required due to the small 

as scales of the structures and the large radial gradients of equilibrium quantities at the pedestal. 

89 Thus, as a consequence of limited computational resources, not all aspects of an experiment can 

90 be described realistically in simulations so far. For instance, simulations with high resolution in 

91 radial and poloidal directions, i.e., with a large number of 2D Bezier finite-elements in the case 

92 of JOREK, render important investigations at more realistic plasma resistivities possible (e.g., 

93 Ref. [8]), but only at a very limited number of toroidal Fourier harmonics. 

94 For this work, a different choice was made: The focus is put on high toroidal resolution. This is 

95 done to investigate the influence of toroidal mode-coupling onto the non-linear evolution of an 

96 ELM. The mode numbers resolved in the simulations are listed in Table 2. All runs resolve the 

97 n = 0, . . . , 16 range but with different periodicities. The relatively high number of toroidal modes 

98 involved limits the possible radial and poloidal resolutions: Most simulations are carried out with 

99 about 5500 Bezier elements. The corresponding finite-element grid is shown in Figure 1. Only 

100 for the simulations with lower plasma resistivity (denoted "eta 6" runs, see next paragraph for 

101 details), the number of Bezier elements is increased by a factor of two. Grid accumulation is 

102 used to increase the resolution radially around the separatrix and poloidally around the X-point. 

103 Due to the comparably low poloidal resolution, only plasma resistivities significantly larger than 

104 in the experiment can be resolved. The respective simulation parameters are listed in Table 3. 
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Figure 1: The flux-surface aligned X-point grid with 5500 Bezier finite-elements is shown. The 
number of grid points are: 96 poloidal points, 40 radial points inside the separatrix, 15 
radial points outside the separatrix, 9 "radial" points in the private flux region, and 9 
grid points along the divertor legs. For the eta 6 simulations, these numbers are all 
increased by a factor of y/2, leading to about 1 1000 Bezier elements. 
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105 The limited poloidal resolution also reduces the growth rate of modes with high mode numbers 

106 artificially. Thus when increasing the poloidal resolution, the most unstable mode number would 

107 shift towards larger n. On the other hand, diamagnetic stabilization is not taken into account in 

108 the simulations. Including this effect would have a stabilizing effect onto high poloidal mode 

109 numbers. The electron diamagnetic frequency for n = 10 is about 10 5 s _1 (calculated at a nor- 

110 malized poloidal flux of = 0.9). This is comparable to the fastest linear growth rates in the 
in simulations (see Section 4). Thus, the most unstable mode numbers would probably be similar 

112 in simulations with higher poloidal resolution and diamagnetic stabilization taken into account. 

113 All simulations are based on typical ASDEX Upgrade discharge parameters, details are given 
in in Section 3.1. The simulations concentrate on the early phase of an ELM-crash up to the 

115 point where filaments start to form. The computations are carried out mostly on the HPC-FF 

116 cluster located in Jiilich, Germany. The eta5 simulations with periodicity 1 and about 5500 

117 Bezier elements require at least 102 compute nodes (8 cores and 24 GB of memory each) due 
us to memory requirements of the solver and take about ten thousand CPU hours to complete. The 
no eta6 computation with 11000 Bezier elements is at the limit of what can be investigated with 
120 JOREK on this machine. 



121 3.1. Physical parameters 

122 A typical ASDEX Upgrade H-mode discharge with type-I ELMs constitutes the basis of the 

123 simulations: Geometry and profiles are taken from discharge 23221 at 4.7 seconds with a 

124 plasma current of IMA, 8 MW of neutral beam injection heating and 1.5 MW of electron 

125 cyclotron resonance heating. The equilibrium reconstruction with the CLISTE code [20, 21] 

126 takes into account measured kinetic profiles. Source terms S p and St and perpendicular diffu- 

127 sivities D± and K± are adjusted such, that the equilibrium does not change significantly with 

128 time. The core temperature is ksTsi = fcjj(7i si + 7i,si) = 12.4 keV. The safety-factor takes a 

129 value of q(0) = 1 in the plasma core and = 0.95) = 4.7 close to the separatrix where 

130 *P;v = OP — axis) / (^separatrix — ^axis) denotes the normalized poloidal flux. A pure deuterium 

131 plasma with a core density of 6- 10 19 m~ 3 is assumed. The heat diffusion anisotropy, Ku/K±, 

132 takes a value of 7 • 10 6 at the separatrix. 

133 The spatial resolution required for the simulation is, amongst others, determined by the resistive 

134 skin depth 8si = y/2r]si/ (HoYsi) which is about 6 mm in eta 6 simulations. As the spatial 

135 resolution possible in the poloidal plane is limited by computational resources, realistic plasma 

136 resistivities with a resistive skin depth of about 0.3 mm cannot be resolved (resolving spatial 

137 scales smaller than the ion gyro-radius is of course not reasonable anyway in MHD-simulations). 

138 The following data are used as inputs for the JOREK simulation: 
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140 



139 



From experimental measurements: Temperature and density profiles, and toroidal mag- 
netic field strength. The pressure profile is shown in Figure 2a. 



144 



141 



142 
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From CLISTE-equilibrium reconstruction based on experimental measurements: FF'- 
profile and the values of the poloidal Flux *F at the JOREK computational boundary. Here, 
F = (2n I Ho) Jp i si is proportional to the poloidal plasma current, 7 po i, and F' = dF / <f *J/. 
The q-profile of the equilibrium is shown in Figure 2b. 



145 In JOREK simulations, the Grad-Shafranov equation is solved first based on these input param- 

146 eters. The equilibrium perfectly agrees with CLISTE (q-profile, flux surfaces, etc.). After that, 

147 an "equilibrium refinement" phase is required where the time-evolution equations are solved 

148 only for the n = mode, with very small time-steps that are gradually increased. This allows 

149 plasma flows to equilibrate [6]. Successively, the reduced MHD equations are evolved in time, 

150 taking into account some or many toroidal Fourier modes depending on the case considered. 

151 Instabilities then develop out of an initially very small random perturbation. 



152 4. Simulation Results 



153 In the following, the simulation results are described and analyzed. Section 4.1 addresses sim- 

154 ulations with low toroidal resolution, while Section 4.2 covers the situation at high toroidal 

155 resolution. In the succeeding Section 4.3, an attempt towards more realistic plasma resistivities 

156 is made. The simulation results are compared to experimental findings in Section 5. 



157 4.1. Low Toroidal Resolution 



158 This section provides simulation results for periodicity 8, where only the toroidal modes n = 0, 8, 

159 and 16 are resolved. A ballooning-like exponentially growing mode located close to the plasma 

160 boundary develops at the low-field side. As seen in the energy diagnostics shown in Figure 3, 

161 the n = 8 mode is linearly more unstable (growth rate /si = 2.0 x 10 5 s _1 ) than the n = 16 

162 mode (/si w 1.5 x 10 5 s _1 ). Due to mode-coupling, the structure of the n = 16 mode changes at 

163 t = 284 \xs in the simulation - the position of its maximum amplitude moves radially from the 

164 q = 4 to the adjacent q = 3.75 resonant surface. Hereby, the growth rate of the n = 16 mode 

165 increases significantly to /si = 4.3 x 10 5 s _1 which is roughly the double n = 8 growth rate. The 

166 n = 8 mode also remains dominant at the onset of non-linear mode saturation (t « 300 (xs). 

167 The ballooning-structure that develops at the whole low-field side of the plasma is shown in 

168 Figure 4 for time point 298 |xs in the simulation. The "density-fingers" are very regular with 
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Figure 3: Time-traces of the magnetic and kinetic energies contained in the individual toroidal 
harmonics are plotted for the simulation with periodicity 8. The n = 8 mode is linearly 
more unstable than the n = 16 mode and also remains dominant when non-linear satu- 
ration sets in. Due to non-linear mode-interaction, the growth rate of the n = 16 mode 
increases significantly at t = 284 |xs. The n = magnetic energy is dominated by the 
toroidal magnetic field which is fixed in time as described in Section 2.1. 
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Figure 4: The density distribution with developing ballooning-structure in the simulation with 
periodicity 8 is shown at 298 (xs. Regular ballooning-structures are observed on the 
whole low-field side. All ballooning-fingers are roughly equally wide in poloidal di- 
rection in the straight-rieldline angle 6* (the white dots indicate equidistant distances 
in d*). 

169 a poloidal size of about 15 — 20 cm at the outer midplane. The poloidal "compression" of the 

170 structures in the vicinity of the lower (active) and the upper (inactive) X-points compared to the 

171 outer midplane is a consequence of field-line stagnation - the poloidal width of the structures is 

172 roughly constant in the straight-fieldline angle 8*. In Figure 4, this can be seen by comparing 

173 the density fingers to the white dots which divide the poloidal circumference into equidistant 

174 sections in 6*. When the exponentially growing perturbation gets visible in the density dis- 

175 tribution, distortions start to build up which propagate into the vacuum region as finger-like 

176 structures with significantly increased density due to the E x B drift. Their radial velocity, mea- 

177 sured by tracing the position at which the density equals 10 percent of the core density, increases 

178 to about 3 km/s and saturates at that level. In the beginning, the density shows sinusoidal ex- 

179 cursions of the density contours which grow over time (linear phase). As the instability grows 
iao and non-linear saturation sets in (energy growth rates start to decrease), the density fingers de- 

181 velop sub-structures. The changing structure also reflects in a different mode-spectrum, where 

182 the n=16 energies get closer to the n=8 energies (Figure 3). The ideal-wall boundary condi- 

183 tions contribute to the saturation of radial velocity when the distance between the mode and the 

184 wall gets significantly smaller than its poloidal wave-length as mirror-currents build up that slow 

185 down the mode-evolution. 
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4.2. High Toroidal Resolution 



187 Now, the same setup as in the previous Section is considered at periodicity 1: All toroidal modes 

188 in the range n = to 16 are resolved. The comparison of simulations with different periodicities 

189 allows to identify effects caused by the coupling between toroidal modes. 

190 Time-traces of the magnetic energies contained in each toroidal harmonic are shown in Figure 5. 

191 Linearly, the n = 10 mode has the largest growth rate /si « 2.0 x 10 5 s _1 . In a similar way as 

192 described for the n = 16 mode in the previous Section, the initially very small growth rate of the 

193 n = 1 mode (ysi ~ 2 x 10 4 s _1 ) suddenly changes at t = 150 \x& due to the non-linear interaction 

194 between the toroidal harmonics and becomes very large: 7si ~ 4 x 10 5 s _1 . In the non-linear 

195 phase of the mode, the n = 1 perturbation reaches a similar magnetic energy as the n = 10 

196 perturbation which remains dominant also at the beginning of non-linear saturation (t « 300 fxs). 

197 A first important effect that cannot be covered in simulations with low toroidal resolution (i.e., 

198 high periodicity) is that low-rc modes can grow to large amplitudes non-linearly. The growth 

199 rate of the dominant mode (n = 10 in our case) is not affected significantly by the toroidal 

200 mode-coupling. Also, the radial propagation velocity of the filaments into the vacuum region 

201 hardly changes compared to the case with low toroidal resolution: The filaments accelerate in 

202 the beginning and saturate at a velocity of about 3 km/s. 

203 The developing density perturbation is shown in Figure 6. Also with high toroidal resolution, 

204 a ballooning-like structure is produced at the low-field side of the plasma. The poloidal size 

205 of the ballooning-frngers is around 10—12 cm at the midplane. In comparison to simulations 

206 with low toroidal resolution, these structures are a bit smaller. A significant difference becomes 

207 obvious when comparing Figures 4 and 6: Due to the mode-coupling, not all fingers grow to 

208 the same amplitude. A cluster of fingers can be seen that develops much stronger than the rest 

209 of the ballooning-structures. A strong localization of perturbations has also been observed in a 

210 ballooning-instability simulated with the BOUT code [22]. 

211 The localization of the perturbation becomes even more obvious when the magnetic footprint 

212 of the mode is considered, in Figure 7, the perturbation of the poloidal magnetic flux is plotted 

213 for simulations with different periodicities. Clearly, the localization of the mode can only be 
2u described correctly in simulations with periodicity 1. Figure 8 shows the perturbation of the 

215 poloidal flux at the outboard midplane versus the toroidal angle. 

216 The perturbation is already localized in the linear phase of the mode. A qualitative change 

217 between the linear and the non-linear phases is shown in Figure 9, where the current perturbation 

218 is plotted for two different time-frames in the simulation with periodicity 1. In the non-linear 

219 phase where the ballooning-fingers become visible in the density perturbation, the previously 

220 alternating current filaments merge at the position of the separatrix around the outer midplane. 

221 Large areas with positive respectively negative currents are created. 
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Figure 5: Time-traces of the magnetic energies contained in the individual toroidal modes are 
shown for the simulation with periodicity 1 . For clarity, kinetic energies are omitted 
and sub-dominant modes are only indicated by dotted gray lines. It is remarkable that 
the n = 1 mode reaches a comparable energy level at the onset of non-linear saturation 
as the n = 10 mode, which is the linearly most unstable mode. 
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Figure 6: The density-perturbation observed in the simulation with periodicity 1 is plotted at 
294 |xs. The ballooning-structures become less regular and perturbations are strong 
only within a localized region. In the cross-section shown, this region is located on the 
upper low-field side. 

222 The strongest perturbations of all physical quantities are localized in a flux-tube like region 

223 which extends from the vicinity of the lower active X-point along magnetic field lines to the 

224 vicinity of the upper inactive X-point (compare Figure 7). The perturbations are strongest around 

225 the midplane. As an exception, vn is perturbed especially around the end-points of this flux-tube, 

226 a consequence of field-line stagnation close to the X-points. However, the radial perturbation 

227 positions differ as shown in Figure 10. It can be seen, that the strongest perturbations of the 

228 magnetic quantities *P and j are located in the region of strong plasma current, while the kinetic 

229 quantities are perturbed further outwards in the region of strong pressure gradients. 



230 4.3. Towards more Realistic Resistivities 

231 This Section briefly shows results for simulations with the plasma viscosity and resistivity re- 

232 duced by a factor of 10 (eta 6 cases) compared to the simulations shown above. To be able to 

233 resolve these more realistic parameters, the number of 2D Bezier elements in the poloidal plane 

234 was increased by a factor of two. These simulations need to be considered with care as the most 

235 unstable mode is n = 13 while we do not take into account mode numbers beyond n = 16 for 

236 computational reasons. 

237 It can be seen that a strong localization of the perturbations is observed at periodicity 1 as 

238 in the eta5-cases. This is shown for the perturbation of the poloidal flux in Figure 11. In 
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Figure 7: Contours of the poloidal flux perturbation are shown for simulations with (a) peri- 
odicity 8, (b) 4, (c) 2, and (d) 1, respectively. The red and blue contours are plot- 
ted at the surfaces corresponding to the perturbed poloidal flux values ^ed/biue = 
±0.7- (|*-Pmin| +*l'max)/2. Here, ^min and *1Vix denote the strongest negative posi- 
tive perturbation values, respectively. At lower periodicities, the perturbation steadily 
gets more localized. 




Figure 8: The perturbation of the poloidal flux at the outboard midplane is shown for the simula- 
tion with periodicity 1 versus the toroidal angle for two transits around the torus. The 
perturbation amplitude shows a strong toroidal variation equivalent to a localization of 
the perturbation to Atp 3 rad (f.w.h.m.). As equilibrium, boundary conditions, and 
sources are completely axi-symmetric, the localization position is essentially arbitrary 
which proves to be true when looking at a set of different simulations. 
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Figure 9: The current perturbation at the plasma edge is shown for the simulation with period- 
icity 1 in the (a) linear (240 \x&) and (b) non-linear phases (300 \x&). In the non-linear 
phase, large regions with positive respectively negative current (this cross-section) 
form at the separatrix (dashed line) around the midplane. 

239 contrast to the eta5-simulations, the perturbation maximum is not located around the midplane 

240 but more towards the top and bottom regions of the low-field side. This distribution of the 

241 flux-perturbation is not an artifact caused by cutting toroidally at n = 16: A simulation with 

242 periodicity 2 was carried out in which the toroidal modes n = 0, 2, . . . , 20, 22 are resolved, where 

243 the strongest perturbation of the poloidal flux is not observed at the midplane but above and 

244 below it, as well. 

245 At ASDEX Upgrade, an off-midplane mode-structure has recently been observed in the tem- 

246 perature using ECE-Imaging [23]. In our simulations, the perturbation maximum of the kinetic 

247 quantities is, however, located around the midplane. This is a consequence of the comparably 

248 large plasma resistivities in our simulations which allow magnetic and kinetic quantities to de- 

249 couple. At smaller resistivities, which we cannot resolve at present, also the kinetic quantities 

250 might show an off-midplane behavior. 



251 5. Comparison to Experiments 



252 In this Section, some properties of the simulated edge instabilities are compared to experimental 

253 observations. This shows that important aspects of the early phase of edge localized modes 
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Figure 10: For the simulation with periodicity 1, the radial positions of the strongest perturba- 
tions are shown at t = 300 \x& for the seven physical variables and are compared to 
profiles of the plasma current and the pressure gradient. 




Figure 1 1 : For simulations with (a) periodicity 2 respectively (b) 1 where the plasma resistivity 
and viscosity is reduced by a factor of 10 compared to the simulations presented 
above, the poloidal flux perturbation is shown analogously to Figure 7. A strong 
localization of perturbations is observed in these simulations as well. 
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t-t E LM [ms] CD MAp [rad] 

Figure 12: (a) Magnetic signals from selected Mirnov-coils are shown for a strongly localized 
solitary magnetic perturbation in ASDEX Upgrade as in Reference [24]. All mea- 
surement locations are mapped to a common toroidal angle 0map via field-line trac- 
ing. The magnetic perturbation propagates with constant toroidal velocity in the elec- 
tron diamagnetic drift direction in the lab-frame as indicated by the red dashed line. 
The onset time of the erosion of pedestal temperature and density profiles is denoted 
*elm- (b) The time-derivative of the magnetic field measured by Mirnov coils is plot- 
ted versus the toroidal mapping angle 0map at t — ?elm = —0.03 ms. The solitary 
magnetic perturbation is localized to A</> « 1.2 rad. 

254 are well described by the reduced MHD model. More detailed comparisons between JOREK 

255 simulations of complete ELM crashes and experimental measurements at ASDEX Upgrade are 

256 planned for the future (e.g., evolution of pedestal gradients, detachment of filaments, heat-flux 

257 patterns at divertor plates). 

258 The poloidal flux perturbation from the simulation with periodicity 1 shown in Figure 8 ex- 

259 hibits a toroidal localization: Large perturbation amplitudes are localized to a region of about 

260 Aip ?a 3 rad. Thus, the modes we observe in our simulations of the early ELM phase when sim- 

261 ulating the full torus (periodicity 1) exhibit a similar magnetic structure as so-called solitary 

262 magnetic perturbations recently discovered at the ELM onset in ASDEX Upgrade and described 

263 in great detail in Reference [24] . From the systematic analysis of a large number of ELM crashes, 

264 a continuous distribution of the mode solitariness was reported between cases with a very pro- 

265 nounced toroidal localization (an example is shown in Figure 12) and cases with a magnetic 

266 perturbation strength that is toroidally virtually uniform. The toroidal localization observed in 

267 our simulations (localized to Atp « 3 rad) is less pronounced than the extreme example of Fig- 

268 ure 1 2b with « 1 .3 rad. A direct comparison is planned for the future making use of a virtual 

269 magnetic diagnostic which determines magnetic signals from the simulations at the same posi- 

270 tions as the Mirnov coils. Toroidally asymmetric structures at ELMs are also described from 

271 experimental observations in References [25-28]. In analytical calculations, localized instabil- 

272 ities were also reported by Wilson et.al. [29]. These "explosive ballooning" instabilities grow 

273 much faster non-linearly than linearly and a poloidal narrowing of the instability in the non- 
274 linear phase is reported. Both features are not observed in the simulated edge instabilities which 
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275 indicates that different mechanisms are responsible for the localization in our simulations. 

276 The dominant toroidal mode-number turns out to be 10 in the simulations. With the plasma 

277 resistivity reduced towards more realistic values, the dominant mode number shifts towards 13. 

278 This is in quite good agreement to experimental findings in the tokamaks ASDEX Upgrade and 

279 MAST for type-I ELMs, where mode-numbers of 8 — 24 were observed in energy deposition 

280 patterns [27], around 15 was found from measurements with the midplane manipulator and 

281 visible-light imaging [30], and mode numbers of 18 ±4 have been obtained for the onset of the 

282 ELM-crash using ECE- Imaging [23]. Uncertainties in our simulations come from the limited 

283 poloidal resolution and the neglect of diamagnetic stabilization as discussed in Section 3. 

284 Low-?? modes gain large amounts of energy non-linearly in our simulations with periodicity 1. 

285 This allows them to interact much more efficiently with core-MHD modes like tearing modes 

286 which typically also feature low toroidal mode numbers like 1 or 2. Indeed, there is experimen- 

287 tal evidence from the DIII-D tokamak that ELMs are an important triggering mechanism for 

288 neoclassical tearing modes [31]. 

289 The poloidal extent of filaments on the outer midplane observed in simulations with high toroidal 

290 resolution is around 10 — 12 cm. Measurements in ASDEX Upgrade and MAST revealed fila- 

291 ment sizes perpendicular to the field lines of 5 — 10 cm [30]. For ASDEX Upgrade, perpendic- 

292 ular and poloidal filament sizes are equivalent due to the small field-line pitch-angle such that 

293 simulation results and experimental observations show good agreement here as well. 

294 In the simulations, the radial velocity of the developing finger structures saturates at about 

295 3km/s after an initial acceleration. This corresponds to a distribution of the radial filament 

296 speed with an upper limit of 3 km/s. This velocity depends on the stability of the initial equilib- 

297 rium. The unrealistically large values for the plasma resistivity might lead to an over-estimation 

298 of the filament speeds, while the ideal-wall boundary conditions tend to reduce the radial veloc- 

299 ity. In experimental measurements, the radial filament speed is found to be distributed around 

300 1 km/s in ASDEX Upgrade [32, 33]. Filament speeds faster than 2 km/s occur in 20% of the 

301 cases in both References and almost no filaments faster than 3 km/s are observed. Hence, radial 

302 filament speeds in simulations and experimental measurements seem to agree reasonably well. 

303 In the magnetic quantities, an off-midplane mode- structure is observed in the simulations with 

304 lower plasma resistivity (eta 6 simulations). As the resistivity is still unrealistically large in 

305 these simulations, magnetic and kinetic quantities are decoupled such that the strongest perturba- 

306 tion of the temperature is located at the midplane. Still, this might be related to the off-midplane 

307 structures observed by ECE-Imaging in ASDEX Upgrade [23]. 
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308 



6. Conclusions and Outlook 



309 Exponentially growing ballooning-like modes have been simulated with the reduced-MHD ver- 

310 sion of the non-linear MHD code JOREK in the geometry and using the profiles of a typical 

311 ASDEX Upgrade H-mode discharge. Dominant toroidal mode numbers, poloidal filament sizes, 

312 and radial filament-propagation speeds of these instabilities are in good agreement with exper- 

313 imental observations for type-I ELMs in ASDEX Upgrade. At sufficient toroidal resolution, 
3H perturbations show a pronounced toroidal and poloidal localization which is compatible with 

315 solitary magnetic perturbations recently discovered in ASDEX Upgrade. In some cases, the 

316 perturbation of the magnetic flux is stronger at the top and bottom low-field side than at the 

317 midplane. Presumably due to a decoupling of magnetic and kinetic quantities caused by the un- 

318 realistically large plasma resistivity, density and temperature perturbations are always localized 

319 on the midplane of the low-field side. Strong perturbations in the \ow-n modes are triggered 

320 non-linearly in the simulations with periodicity 1 and might explain the strong interaction of 

321 ELMs with core-MHD modes like neoclassical tearing modes observed in some experiments. 

322 While this work concentrates on the early phase of an ELM, further studies are planned to com- 

323 pare the simulation of a full ELM crash to experimental observations requiring a more sophisti- 

324 cated modeling of the scrape-off layer. Simulations of a full ELM cycle will also be attempted. 

325 Future numerical improvements and increased computational resources will be used to advance 

326 our investigations towards more realistic plasma parameters while keeping high toroidal resolu- 

327 tions. 
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334 A. Derivation of the Induction Equation 



The reduced MHD equations implemented in the JOREK code can be derived following Refer- 
ence [16]. For the induction equation, this is shown in the following. The starting points are the 
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well-known expression for the electric field in the MHD approximation, 

E = -vxB + Tjj (10) 

and the Maxwell-Faraday law expressed in the vector potential A, 
dA 

= -E-F Vu. (11) 

at 

Inserting (10) into (11) and multiplying it with the toroidal unity vector t>0 yields 

— =Tj/+(vxB)-e^-F Vi < -e^, (12) 

where the poloidal flux is identified as the major radius times the toroidal component of the 
vector potential, *P = R A • e^, and j = — j ■ denotes the toroidal plasma current. Using Equa- 
tions (8) and (9), this can be written as 

du 

- = ni -R [u ,V\-F 0w (13) 

335 which is the induction equation (Equation (3)) solved in the JOREK reduced MHD model with 

336 the Poisson bracket [a, \P] = f| |^ — || ^ . In the last step, the reduced MHD approximation to 

337 first order in e = V^/V^ <C 1 yielding v po i « Vj_ was applied. 



338 The poloidal components of Equation (11), obtained by applying the operator x to this equa- 

339 tion, yield a definition equation for the poloidal velocity (see poloidal components of Equa- 

340 tion (8)) in which u can be identified as the poloidal velocity stream function. In this set of 

341 equations, u also acts as electric potential (except for a constant factor Fq). 



342 Galilei-invariance of the induction equation (Equation ( 1 1)) is not obvious at first glance. How- 

343 ever, the proof is straightforward when taking into account that the scalar potential = Fqu is 

344 modified according to <p — > <p — vo • A under the transformation v — > v — vo while the vector po- 

345 tential remains unchanged (non-relativistic limit). In the large aspect-ratio limit, it can also be 

346 shown easily that the reduced-MHD induction equation (Equation (13)) is invariant to a trans- 

347 formation v — > v — \ z with z along the cylinder axis, as the scalar potential transforms according 

348 tO (f) — > — Vj*P. 
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